Ab initio simulations of α- and β-ammonium carbamate (NH4·NH2CO2), and the thermal expansivity of deuterated α-ammonium carbamate from 4.2 to 180 K by neutron powder diffraction

The results of ab initio simulations carried out using density functional theory on two polymorphs of ammonium carbamate are reported. Also reported is a refined crystal structure of deuterated α-ammonium carbamate from high-resolution neutron powder diffraction data, with thermal expansion measurements from 4.2 to 180 K.


Introduction
With recent robotic missions to icy planetary bodies in our solar system, most notably the flybys of Pluto, Charon and Arrokoth by the New Horizons spacecraft, and of Ceres by the Dawn spacecraft, the importance of characterizing 'ices' composed of primordial mixtures of ammonia, carbon dioxide and water has become more urgent. Although long suspected to be present in the outer solar system (Kargel, 1992), solids in the ternary NH 3 + CO 2 AE H 2 O system ( Fig. 1), which includes two phases of ammonium carbamate (NH 4 ÁNH 2 CO 2 ), ammonium carbonate monohydrate [(NH 4 ) 2 CO 3 ÁH 2 O], ammonium sesquicarbonate monohydrate {(NH 4 ) 4 [H 2 (CO 3 ) 3 ]ÁH 2 O}, ammonium bicarbonate [(NH 4 ) 2 HCO 3 ] and the water-deficient compound urea [CO(NH 2 ) 2 ], appear increasingly likely to occur as 'rock-forming' minerals.
Pluto has extensive deposits of N 2 , CH 4 and CO ices across its surface in various mixtures (Protopapa et al., 2017;Schmitt et al., 2017;Moore et al., 2016) and a thin N 2 -dominated atmosphere (Gladstone & Young, 2019), with both Pluto and its moon Charon manifesting the spectroscopic signature of ammoniated material (Dalle Ore et al., 2018Grundy et al., 2016). The dwarf planet Ceres also manifests spectroscopic signatures of ammoniated species, in the form of phyllosilicates (Ammannito et al., 2016;Raponi et al., 2019), which are seen in the regions with high albedo, in particular in the complex crater Occator (De Sanctis et al., 2015, suggesting an NH 4 -rich interior. The Kuiper belt object Arrokoth has abundant CH 3 OH on the surface , with other potential ices, including H 2 O, NH 3 , H 2 CO and CO 2 , in trace amounts (Grundy et al., 2020). In other words, all the constituents necessary to make ammonium carbonates (and related substances) are present and so are likely to form these materials wherever partial melts of ammonia-water are exposed to solid or gaseous CO 2 .
In the same context, ammonium carbonates may have some astrobiological relevance since removal of water leads -via ammonium carbamate -to urea. It is known from laboratory analogue experiments that the organic molecules produced photochemically in the atmosphere of Saturn's giant satellite Titan may be hydrolysed in aqueous ammonia to form both urea and amino acids (Poch et al., 2012). Further hydrolysis of urea would be expected to form ammonium carbonates; on Earth, this process is mediated in soils with the aid of bacterial urease, whereafter the carbonate breaks down to ammonia and water. Meteorite impacts into Titan's surface would provide the requisite liquid (as impact melt) to hydrolyse any solid organics (Artemieva & Lunine, 2003) and both urea and ammonium carbamate/carbonate may be substantial byproducts, persisting on geological timescales at the surface temperature of 95 K.
Not only are ammonium carbonates relevant to the outer solar system, but they also play a role in industry here on Earth; ammonia-water solutions may be used to capture CO 2 from flue gases, deposited in solid form as ammonium carbonates (Han et al., 2013). Furthermore, there has been some interest in using ammonium carbonates for hydrogen storage in fuel cells (Lan et al., 2012).
There has nevertheless been comparatively little interest in studying the structures and properties of these compounds; indeed, the structure of ammonium carbonate monohydrate was determined comparatively recently (Fortes et al., 2014), and the structure of the sesquicarbonate was determined in 2003 (Margraf et al., 2003). Furthermore, the identity of a second () phase of ammonium carbamate was first reported in 2007 (Kuhn et al., 2007), while the only published crystal structure of -ammonium carbamate was that determined in 1973 by single-crystal X-ray diffraction (Adams & Small, 1973). Both ammonium carbamate phases are orthorhombic: -ammonium carbamate has the space group Pbca, with Z = 8, while -ammonium carbamate has the space group Ibam, with Z = 8.
It is surprising that ammonium carbamate has apparently been of little interest, given the similarity of the carbamate ion to the primary amide group; previous work on ammonium carbamate has mainly focused on the thermodynamics and kinetics, typically related to the formation of urea, rather than crystallographic studies. Whilst urea has been extensively studied, the same is not true of the carbonate/carbamate family of compounds.
The goals of this study were to calculate the bulk properties of both phases of ammonium carbamate as a function of pressure, preparatory to the development of a planetary model. Without knowing how something as fundamental as the density varies with pressure and temperature, as well as high-pressure polymorphism, the development of accurate structural and evolutionary models of icy planetary interiors is severely limited. To determine the necessary physical properties includes calculating the bulk elastic properties as a function of pressure using density functional theory (DFT); in particular, how the different arrangement of structural units and hydrogen bonds affects the bulk modulus of the two polymorphs.
An opportunity to test the anisotropy of the -polymorph arose during an attempt to measure the thermal expansion of ammonium carbonate monohydrate. During a synthesis of that product, ammonium carbamate was formed instead, and neutron powder diffraction data were collected as a function of temperature. The thermal expansion derived from these data provides a useful measurement of the structural anisotropy for comparison with the calculations.
Here, therefore, we present the first ab initio calculations to determine the high-pressure behaviour of bothand -ammonium carbamate; the first neutron powder diffraction data of -ammonium carbamate, thereby providing the first accurate H-atom positions to allow characterization of the hydrogen bonding in the crystal; and thermal expansivity measurements from 4.2 to 180 K. These results begin to bridge the current gap in our knowledge of the behaviour of ammonium carbonates under conditions related to the icy moons of the outer solar system.

Experimental
Time-of-flight neutron powder diffraction data were collected using the 'High Resolution Powder Diffractometer' (HRPD) at the STFC ISIS neutron spallation source of the Rutherford-Appleton Laboratory, UK. The HRPD instrument was chosen due to its high resolution (Ád/d ' 8 Â 10 À4 in the backscattering bank), allowing the highest precision in refined unitcell parameters and accurate partitioning of intensity in wellresolved high-Q Bragg peaks.

Sample preparation
Ammonium carbamate was formed serendipitously in the process of attempting to synthesize deuterated ammonium carbonate monohydrate, during which a 28 wt% ND 3 solution was exposed to CO 2 from sublimating dry ice pellets, forming a slurry of white material. Consequently, the specimen contained excess water in the form of water ice, which also acted as a reservoir for small amounts of excess ammonium. The resulting white compound was stored in a freezer. The sample was taken out of the freezer and allowed to warm and soften, after which the partially molten slurry was transferred into a cryomortar under liquid nitrogen and ground to a fine powder using a similarly cooled pestle. Once fully ground, the material was loaded into an aluminium-framed slab-geometry can of inner dimensions 23 Â 18 Â 15 mm (h Â w Â d) relative to the incident neutron beam, also precooled in liquid N 2 . Vanadium foil windows were indium-sealed to the front and back faces of the slab can, with the exposed components on the incident-beam side being masked with Gd foil. The sample holder was wired with an RhFe resistance thermometer and a 30 mm cartridge heater, after which, it was immersed in liquid nitrogen for transfer into a He cryostat operating at 100 K. Once oriented in the beam, the cryostat temperature was reduced to 4.2 K for measurements to commence.

Data acquisition
Time-of-flight data on HRPD are measured routinely in a range of 100 ms bandwidth 'windows' where disc choppers act as bandwidth selectors, these typically being 30-130 or 100-200 ms. In the instrument's backscattering detector bank (2 = 154-176 ), these time-of-flight ranges provide access to dspacings in the ranges 0.65-2.60 (30-130 ms) and 2.15-4.00 Å (100-200 ms) at the highest resolution. Intermediate-resolution data are simultaneously measured in the 90 detector banks (0.85-3.90 Å in 30-130 ms and 2.82-6.00 Å in 100-200 ms). Since HRPD has an essentially Q-independent resolution function, we can obtain high-precision unit-cell parameters from data measured only in the 100-200 ms window whilst avoiding peak overlap from the additional peaks at shorter d-spacings. For accurate structure refinement, however, it is more desirable to measure the shorter d-spacing peaks in the 30-130 ms window.
Consequently, we collected a low-noise structural data set at 4.2 K in the 30-130 ms window, counting for 6.5 h (corresponding to an integrated proton beam current of 250 mA h) and then measured on warming using only the 100-200 ms window. These thermal expansion data were collected in 5 K increments from 10-150 K and (due to time constraints) in 10 K increments up to 180 K. Counting times were 15 min, corresponding to an integrated proton beam current of 15 mA h per step. The diffraction data were focussed, normalized to the incident beam spectrum and corrected for detector efficiency by reference to a vanadium standard using Mantid (Arnold et al., 2014) and then exported in a format suitable for analysis with GSAS/Expgui (Larson & Von Dreele, 2004;Toby, 2001).

Structural refinement
In the 4.2 K neutron powder diffraction dataset, two phases were identified by inspection of the peak positions -deuterated water ice and -ammonium carbamate. The Rietveld method was used to refine the structure of D 2 O ice Ih and -ammonium carbamate, starting from the previously published atomic coordinates based upon the single-crystal X-ray investigation of Adams & Small (1973). In this refinement, the unit-cell parameters, scale factor, peak-profile coefficients (GSAS profile function #3) and a twelve-term Chebyshev polynomial background were allowed to vary and the atomic coordinates and isotropic displacement parameters (U iso ) were refined independently, without restraints. It was found necessary to refine a wavelength-dependent absorption/ reflectivity correction (GSAS model #1). The final fit to the data is excellent, the R wp value being better than 2%; a Neutron powder diffraction pattern acquired at 4.2 K in the highest resolution backscattering banks of HRPD. Red symbols represent the measured data, green lines the result of Rietveld refinement, pink lines the difference profile and the tick marks the expected positions of each Bragg reflection for ND 4 ÁND 2 CO 2 (black) and D 2 O ice Ih (red). graphical depiction of the fit to the neutron powder data is given in Fig. 2. The refined structural parameters are given in the crystallographic information file (CIF) in the supporting information.
At each temperature step, the unit-cell parameters of both phases were determined using the Le Bail method (with initial intensities based on the GSAS F calc method); in addition to the unit-cell parameters, the peak profile parameters and the background parameters were allowed to vary. Consideration of the axial ratio of the D 2 O ice Ih present in the sample (see Fig. 3 and x4.1) indicated that, at 180 K, it was sufficiently pure to allow its use as an internal standard. To ensure uniformity with previous work, the unit-cell parameters of D 2 O at 180 K were therefore fixed to those measured in a previous experiment with a silicon standard (Fortes, 2018), and the diffractometer constants subsequently refined.

Computational methods
To provide further information for planetary models, such as high-pressure behaviour, and to provide an accurate structure for -ammonium carbamate, both structures of ammonium carbamate were input into CASTEP (Clark et al., 2005). Static density functional theory (DFT) calculations were then carried out using the PBE pseudopotential (Perdew et al., 1996) with the Tkatchenko-Scheffler dispersion correction (Tkatchenko & Scheffler, 2009). The importance of dispersion corrections for the study of elastic properties and anisotropy in planetary ices was demonstrated recently by Meusburger et al. (2021). For the work reported here, we had found that geometry optimizations completed without dispersion corrections led to disagreements of $10% in the unit-cell parameters, these being reduced to 1-3% with the use of dispersion corrections. Similarly, the estimated phase-transition pressures reported later were reduced from $2 to $0.4 GPa.
Convergence tests were carried out to optimize the k-point sampling of the Brillouin zone within the Monkhorst-Pack scheme (Monkhorst & Pack, 1976) and the kinetic energy cutoff of the plane-wave basis set. For -ammonium carbamate, a converged grid of 2 Â 6 Â 5 ($0.03 Å À1 reciprocal lattice spacing), with an energy cut-off of 1000 eV, yielded a total energy convergence better than 10 À3 eV. For -ammonium carbamate, a k-point mesh of 2 Â 2 Â 3 ($0.05 Å À1 reciprocal   lattice spacing), with an energy cut-off of 980 eV, yielded a total energy convergence better than 10 À3 eV.
A series of calculations were then performed in which the ions were allowed to move according to the Hellman-Feynman forces with the symmetry maintained, and the unitcell shape was allowed to vary; the initial structural parameters of -ammonium carbamate were those obtained from our neutron powder diffraction data at 4.2 K and of -ammonium carbamate were taken from Kuhn et al. (2007).

Behaviour of D 2 O ice Ih upon warming
When comparing D 2 O ice Ih unit-cell parameters to literature values, there is a significant difference in the ice unitcell parameters compared with pure ice at low temperatures ( Fig. 3). Although the unit-cell volumes at 180 K agree well, there is a growing divergence at lower temperatures, such that the unit-cell volume is smaller than that of pure D 2 O ice Ih at 10 K by 0.11% and the transition to negative volumetric thermal expansion is shifted to a lower temperature. Furthermore, the deviation is associated almost entirely with the behaviour of the crystal's a axis; this is most apparent from a plot of the c/a ratio against temperature (Fig. 3). From 180 to 135 K, the c/a ratio is in excellent agreement with both the absolute values and the temperature dependence found in pure D 2 O (Fortes, 2018); below 135 K, the c/a ratio increases systematically.
The dependence of the c/a ratio in ice Ih on both thermal history and chemical doping was described by Fortes (2019), who attributed the observed effects to changes in partial order in the crystal. We therefore propose that small concentrations of ammonium ions present in our D 2 O ice Ih sample effects similar very subtle structural changes that drive the c/a ratio to adopt a nearly ideal hexagonal close-packed value (1.630).
In addition to the likely very small quantities of ammonium ions trapped in the crystal structure of ice Ih, it appears that Neutron powder diffraction patterns acquired at 160 K (top) and 170 K (bottom) in the highest resolution backscattering banks of HRPD. Red symbols represent the measured data, green lines the result of Rietveld refinement, pink lines the difference profile and the tick marks the expected positions of each Bragg reflection for: (a) ND 4 ÁND 2 CO 2 (black); (b) D 2 O ice Ih (red); (c) ND 3 ÁD 2 O (blue); (d) ND 3 Á2D 2 O (green). At 160 K, peaks of ND 3 ÁD 2 O appear, and at 170 K, peaks of ND 3 Á2D 2 O also appear; by 180 K, these peaks have disappeared. Table 1 Comparison of bond lengths, polyhedral volumes, various distortion metrics and hydrogen-bond parameters (Å , ) in -ammonium carbamate.
The values are from the outputs of CASTEP and measurements on HRPD, a powder X-ray diffraction experiment (Baisch et al., 2006) and a single-crystal X-ray diffraction experiment (Adams & Small, 1973). Note that the HRPD sample was deuterated. Atomic labels in Baisch et al. (2006) correspond to: H11 = H3, H12 = H2, H22 = H5, H23 = H6, H24 = H4 and N2 = N1, according to the labelling of Adams & Small (1973). The distortion index is based on bond lengths in the polyhedra and was defined by Baur (1974). 'Quadratic elongation' gives a quantative measure of polyhedral distortion, which is independent of the effective polyhedron size (Robinson et al., 1971), and 'bond-angle variance' gives a quantative measure of how the bond angles in the tetrahedra deviate from the ideal value of 109.28 -lower values equate to bond angles closer to the ideal (Robinson et al., 1971 (7) 158 (5) there was also sufficient excess ammonia present -likely as an amorphous aqueous solid in pores and on grain boundariesto form ammonia hydrates. These were seen in the diffraction patterns during warming of the sample, with the first impurity phase appearing at 160 K. Peaks from this phase became larger at 170 K and were joined by peaks from another accessory phase (Fig. 4). During the final measurement at 180 K, all impurity peaks had disappeared completely. We matched the peaks from both accessory phases; the first to appear was ammonia monohydrate (ND 3 ÁD 2 O) and the second was ammonia dihydrate (ND 3 Á2D 2 O). Since these are both known to melt between 170 and 180 K, it is not surprising that they are absent from the 180 K diffraction data.
The presence of excess water in the form of ice meant we were obliged to keep the sample below the melting point of water. The addition of even small amounts of water from moisture in the atmosphere to ammonium carbamate is enough to slowly form ammonium bicarbonate. With the presence of ammonium ions in the sample, eventually melting and forming ammonium hydrates, the melting point of water would be significantly reduced, and we could not be certain as to when the sample would start to decompose, so it was decided to halt the experiment at 180 K to be certain the sample did not melt. Hydrogen bonding between centrosymmetric pairs of carbamate ions creates an R 2 2 (8) ring motif, measuring 4.286 Å , between centrosymmetric pairs of C atoms; these rings are typically seen in primary amide structures, with amide molecules, R-CONH 2 , forming the rings (Fig. 6).

a-Ammonium
There are many different possible arrangements of these rings within a crystal structure, including various ring, tape, and sheet-like motifs (Leiserowitz & Hagler, 1983). In ammonium carbamate, layers of these hexagonal rings are alternately arranged perpendicular to one another, known as a herringbone motif arrangement. Indeed, the space group of ammonium carbamate, Pbca, can generate an interlayer herringbone motif between the terminal atoms of the amide residue groups, R, which lie on the outside of the bilayer, due to the b-glide plane (Leiserowitz & Hagler, 1983). The form of the ammonium carbamate motifs in -ammonium carbamate is therefore similar to the shallow glide motif in primary amide structures.
In the shallow glide motif, the spacing between layers of the same orientation is typically 3.4 Å ; in ammonium carbamate, this spacing is almost doubled in length, at 6.436 or 6.681 Å (equivalent to the b and c axes, respectively), due to the interlayer motif acting as an additional layer. Three C atoms of  The R 2 2 (8) carbamate ionic dimers seen in the structure of -ammonium carbamate (ND 4 ÁND 2 CO 2 ). Faded tetrahedra are below the plane of the centrosymmetric ring structure and the dashed lines represent H(D) bonds.

Figure 5
The structure of -ammonium carbamate, viewed parallel to the b axis. The grey dashed lines are H(D)Á Á ÁO hydrogen bonds, and the solid black line is the unit cell. separate hexagonal ring structures usually have a separation of $10 Å , but here the distance is 9.27 Å ; this is accommodated by the tilt angle of the ring structures being $68 rather than a typical 61 (Leiserowitz & Hagler, 1983).
In terms of bonds and bond lengths (see Table 1), the length of the N2-D1 bond is longer than the N2-D2 bond, due to the D1 atom donating a hydrogen bond [forming part of the R 2 2 (8) H(D)-bonded motif]. The N2-D2 bond was suggested not to be involved in a hydrogen bond by Adams & Small (1973); it appears that this H(D) atom may be sharing an interaction with the O1 atom, although the (H)D2Á Á ÁO1 distance is longer and more distorted than the remaining H(D) bonds -the D2Á Á ÁO1 bond length measures 2.405 (5) Å , with an angle of 150.4 (4) . This work shows that the D1Á Á ÁO1 and D6Á Á ÁO1 hydrogen bonds are slightly longer than the remaining two bonds, agreeing with the work of Adams & Small (1973).
4.2.2. Thermal expansion of a-ammonium carbamate. The thermal expansion was determined for deuterated -ammonium carbamate from the unit-cell parameters collected on HRPD at 4.2 K, in 5 K increments from 10 to 150 K, and in 10 K increments from 150 to 180 K.
In order to model the thermal expansion to obtain useful thermoelastic parameters, a second-order Grü neisen-Debye model, with a single Debye temperature [equation (S2) in the supporting information], was used in least-squares fitting to the unit-cell volume, with the values V 0 , D , Q and b freely refined (further information regarding Grü neisen-Debye models can be found in the supporting information). The results of this fit are given in Table 2, and shown in Fig. 7. The calculated thermal expansion coefficients are shown in Fig. 8.
The thermal expansion appears to be dominated by lowfrequency vibrational modes with characteristic temperatures in the range 200-230 K; this suggests that the phonon density of states has a cut-off at $17.  Table 2 Parameters obtained by fitting a second-order Grü neisen-Debye model with a single Debye temperature [equation (S2)] to the unit-cell volume of -ammonium carbamate between 4.2 and 180 K.

Figure 7
Unit-cell parameters (a)-(c) and unit-cell volume (d) of deuterated -ammonium carbamate from 4.2 to 180 K. The solid symbols are from Le Bail refinements, with the solid red lines indicating the exponential model [equations (1) and (2) wavenumbers $140-160 cm À1 or $4.2-4.8 THz, respectively). The fit of the Grü neisen-Debye model with a single Debye temperature also suggests that the heat capacity of -ammonium carbamate does not approach the classical Dulong and Petit high-temperature limit (3nk B ), but instead tends towards a high-temperature limit of $2nk B . The value of K 0 / found from this Grü neisen-Debye fit to the unit-cell volume is 49 (1) GPa. This value is much greater than that calculated from fitting an equation of state to the outputs of the DFT calculations (see x4.2.4). The reason for this is most likely that a Debye model with a single Debye temperature does not properly describe the internal energy of a molecular solid such as -ammonium carbamate (see supporting information). An alternative explanation is that 6 ¼ 1; if a value of K 0 is assumed from the high-pressure DFT simulations, where K 0 = 16.5 GPa, equation (S2) gives % 0.33 (1). It could also be that -ammonium carbamate shows a temperature dependence of , which is assumed to be constant in this approximation.
Although Grü neisen-Debye models can provide useful thermoelastic parameters, they can be difficult to evaluate (due to their complexity and, in this case, lack of heat capacity data) and are only dimensionally correct for the unit-cell volume. Therefore, the thermal expansion data were also determined by fitting the unit-cell parameters using a linear combination of a power law and an exponential function: where p and r are scaling terms, and q and s are a 'characteristic temperature' indicative of the saturation behaviour of the function. The temperature dependence of the unit-cell parameters, or unit-cell volume, is then given by: where X(T) is either an axial length or the unit-cell volume. The results (Table 3) have small uncertainties; furthermore, they lack the oscillatory behaviour of polynomial functions and give the correct behaviour at temperatures close to 0 K. As can be seen in Fig. 7, despite its ad hoc nature, the exponential model does an excellent job of fitting the experimental results, showing only extremely small differences from the Grü neisen-Debye model (nevertheless, we present the fitting of different Debye models in the supporting information).  The linear expansivities and the volumetric expansivity, derived from equation (1), are shown in Fig. 8. The symbols are obtained as point-by-point derivatives of the experimental unit-cell parameters, providing a useful visual indicator of the uncertainties due to their scatter.
There is a large degree of anisotropy in the structure upon warming; the greatest expansion is along the b axis, with the a and c axes having similar expansivities. The values of the volumetric thermal expansion coefficient are comparable with water ice and ammonium carbonate monohydrate at the same research papers Acta Cryst. temperature (Fortes, 2018;Fortes et al., 2014). We find no evidence for negative thermal expansion in the structure.
For making a simple density calculation as a function of temperature for planetary modelling, the density can be described using only the exponential function of equation (1) (as was sufficient to describe the volumetric expansion) and using: with the same values of r and s as determined from the fit to the unit-cell volume, and 0 = 1523.87 (2) kg m À3 .

DFT simulations of the structure of a-ammonium
carbamate. The unit-cell volume at ambient pressure and 4.2 K, as measured on HRPD, is 733.126 (9) Å 3 . The zeropressure athermal relaxation of -ammonium carbamate from CASTEP yielded a unit-cell volume of 743.772 Å 3 , which equates to an approximately 1.5% overestimation of the measured value (Table 4). Although the simulated -ammonium carbamate is protiated, the differences from the experimental structure, which is deuterated, are likely to be minor. This difference in unit-cell volume can largely be attributed to the angle between the two herringbone-motif carbamate ions -the simulated value has this close to 90 , whereas the refined structure at 4.2 K has an angle greater than 90 . This means that the b axis has become longer and the c axis shorter in the simulated unit-cell parameters, resulting in the over-and underestimation of the individual unit-cell parameters. The result of relaxing -ammonium carbamate made no significant changes to the structure. The bond lengths in the structure are only marginally different between the simulated and measured structures, i.e. the length of the carbamate ion is very similar, and the distance between centrosymmetric pairs of ions is 4.254 Å for the simulated and 4.304 (5) Å for the measured structures.
One thing to note is the distortion index of the ammonium tetrahedra in the HRPD sample measured at 4.2 K, which is much larger than the other values in Table 1, but remains an order of magnitude smaller than that of Adams & Small (1973). Apart from the distortion of the tetrahedra, the work on HRPD and the athermal zero-pressure DFT relaxation shows a very good agreement with the bonding geometry, such as the H(D)1Á Á ÁO1 and H(D)6Á Á ÁO1 hydrogen bonds being slightly longer than the remaining three bonds [excluding H(D)2Á Á ÁO1], and also agrees well with the work of Adams & Small (1973).

Simulations of the high-pressure behaviour of a-am-
monium carbamate. To provide information relevant to simulating -ammonium carbamate in models of icy planetary interiors, the high-pressure behaviour was calculated by a series of relaxations over a pressure range between À0.9981 and +14.9952 GPa.
To quantify elastic parameters, an equation of state was fitted to the total energy against unit-cell volume [E(V)] and pressure against unit-cell volume [P(V)] curves (BMEOS3; Birch, 1952). These plots are given in Fig. 9, with the parameters from fitting to the calculated unit-cell parameters of -ammonium carbamate given in Table 5. The values of X 0 , K 0 , K 0 0 and E 0 were allowed to freely refine. A parametric form of the Birch-Murnaghan equation of state was used to fit the compressibility along each of the unit-cell axes (for b 3 and c 3 ; the a axis is not fitted, which is explained next).
For making a simple density calculation as a function of pressure, the density was fitted using the equation of state used to fit P(V) points, where X = / 0 (see supporting information). The fitted parameters are: 0 = 1394.2 kg m À3 (value fixed), K 0 = 15.4 (3) GPa and K 0 0 = 5.2 (1). There is a high degree of anisotropy in the structure when compressed and, although not immediately apparent from examination of the E(V) and P(V) curves, there are some interesting and complex features in how the structure deforms in response to an external pressure, in particular along the a axis.
The initial loading shows subtle changes to the carbamate ions, manifesting in changes mostly along the a axis up to pressures of 5 GPa. During the initial loading, the carbamate ions reconfigure, particularly one of the H atoms on the N2 amino group. The N2-H1 and N2-H2 bonds shorten and lengthen, respectively, but in a rather smooth fashion, with the angles C1-N2-H1 and C1-N2-H2 also changing smoothly. This results in the carbamate ions themselves becoming more planar up to a pressure of 2 GPa before transitioning to a decrease in planarity (Fig. 10). The C1-O1  Table 4 A comparison of the unit-cell parameters of -ammonium carbamate from the DFT athermal zero-pressure relaxation using CASTEP, from an X-ray single-crystal diffraction (SXD) experiment at room temperature (Adams & Small, 1973) and from a powder X-ray diffraction experiment (PXRD) at room temperature (Baisch et al., 2006).  Parameters obtained by fitting of equations (1) and (2) to the unit-cell parameters and unit-cell volume of -ammonium carbamate between 4.2 and 180 K.
a axis 17.0517 (2) --0.268 (8) À51 (2) 0.99827 b axis 6.43539 (9) 0.95 (1) À14.5 (3) --0.99973 c axis 6.68077 (6) --0.74 (2) À113 (3) 0.99903 Volume (Å 3 ) 733.148 (9) --1.65 (1) À60.9 (6) 0.99993 bond lengthens slightly, before shortening again. After this initial reconfiguration, the changes are mainly due to the HÁ Á ÁO hydrogen bonds in the rings changing their behaviour, initially lengthening before shortening. Interestingly, these changes only manifest themselves in the lengths of the a and b axes, with barely any effect on the c axis. This compressive behaviour continues up to a pressure of 5 GPa. However, along the a axis, between 5 and 8 GPa (region II), there is a small degree of negative linear compressibility (NLC), followed by an inversion to normal positive linear compressibility (PLC), averaging over this pressure range a small degree contraction along the a axis. The cause of this behaviour is a small rotation of the NH 4 + tetrahedra to accommodate the closing of the carbamate rings, causing the carbamate ions to move outwards along the a axis.
The largest NLC region in the a axis, however, occurs over the range 8-11 GPa (region III), continuing to a lesser extent up to $15 GPa (region IV), but this does not seem to translate into a significant change in the unit-cell volume, as can be observed in Fig. 10. The NLC apparent in the a axis manifests itself in the b and c axes in different ways: between 9.0 and 10.5 GPa, there is a marked change in slope of the b axis, which occurs at the same point as the end of the major part of the NLC; however, along the c axis, there is no obvious change in slope.
The cause of this large NLC is once again from the rotation of the NH 4 + tetrahedra, forcing the carbamate ions to move apart along the a axis. Alternating chains of ammonium tetrahedra and carbamate ions along b likely reach a critical point, in that they cannot move closer together along either  the a or the b axis. The only way to accommodate the increasing stress is therefore to rotate the ammonium tetrahedra. As the carbamate ions had rotated at a lower pressure and had been packed closer together by the contraction along b, one of the O atoms was moved closer to the ammonium tetrahedra, leading to the formation of a weak hydrogen bond. This is evident from the transfer of the HÁ Á ÁO hydrogen bonds from the ammonium tetrahedra to the carbamate ion (Fig. 10). As the NLC begins, the N1-H6Á Á ÁO1 hydrogen bond becomes highly strained because of the ammonium tetrahedra rotating, and instead the bond changes to the O2 atom of the carbamate molecule. This also causes the change in slope of the b axis.
The final transition at a critical pressure of 11.25 GPa (regions III and IV) is caused by the system reconfiguring to the rotated ammonium tetrahedra, with the carbamate ions forced to translate along the a axis.
4.3. b-Ammonium carbamate 4.3.1. DFT structure of b-ammonium carbamate. The zeropressure athermal geometry optimization of the -phase from CASTEP gave a unit-cell volume of 711.6 (2) Å 3 ; although this equates to an approximately 3% overestimation of the measured value at 298 K (Kuhn et al., 2007) (Table 6), the overestimation of the unit-cell volume is likely much larger, since, intuitively, it is expected the structure will contract upon cooling. An experiment at low temperature is needed to confirm the degree of overestimation of the unit-cell volume. However, since no experimental data are currently available, we present comparisons between the athermal simulation and the room-temperature experiment, where the differences can largely be attributed to the bond lengths; DFT calculated bond lengths for -ammonium carbamate are compared to the experimental values in Table 7.
The N-H bond lengths given in Kuhn et al. (2007) are typically short due to the use of X-ray diffraction. The results of the present work, therefore, improve on the accuracy on the H-atom coordinates and are thus a much better determination of the hydrogen-bond lengths in the structure. Once again, this allows for more definitive statements regarding the structure and bonding to be made.
The previously published structure of -ammonium carbamate indicated that it does not form the centrosymmetric dimer ring motifs R 2 2 (8) seen in -ammonium carbamate; instead, it forms one-dimensional chains along the mirror plane by the hydrogen bonding of a nitrile H atom to one of the carbonyl O atoms (H4Á Á ÁO2) at c/2 intervals. The remaining nitrile H atom (H3), in this instance, is unlikely to form a hydrogen bond, due to the N2-H3Á Á ÁO2 hydrogen bond having a small angle (115 ) and long length ($3 Å ). Each carbonyl O atom has two out-of-plane hydrogen bonds to two separate ammonium tetrahedra, which lie either side of the mirror plane, with the tetrahedra midway between the chains.
However, in the DFT-relaxed structure (Fig. 11), the carbamate molecules have rotated, resulting in the H atoms changing position. The N2-H3Á Á ÁO2 hydrogen bond has become quite linear, with an angle of $172 , and is much shorter, at 2.074 Å . The relaxing of the H atoms has altered the bonding geometry so that, now, -ammonium carbamate forms the same hexagonal ring motif, a double two-centre hydrogen-bonded R 2 2 (8) ring that is seen in the -phase, albeit with a different motif layout to the latter (see Fig. 12); the motifs instead lay on mirror planes, but are offset by the two glide planes in the structure. These planar R 2 2 (8) rings are joined via ammonium ions at z = 1 4 and z = 3 4 . It can also be seen in the experimental values that the ammonium tetrahedra are extremely distorted (  Fits to equation (S5) were achieved using only those pressure values that were greater than 0 GPa. For the b and c axes, the values fitted were b 3 and c 3 . No equation of state could be fitted for the a axis due to the presence of negative linear compression.

4.3.2.
Simulations of the high-pressure behaviour of b-ammonium carbamate. As with -ammonium carbamate, the high-pressure behaviour of -ammonium carbamate was simulated to provide information relevant for planetary models. A series of static relaxations were performed over a pressure range between À1.9977 and +10.025 GPa.
To gain information on the bulk elastic parameters, E(V) and P(V) curves were fitted with a BMEOS3 equation of state, with the plots given in Fig. 13 and values of the fit in Table 8. The values of X 0 , K 0 , K 0 0 and E 0 were allowed to freely refine. Plots of the unit-cell parameters as a function of pressure are given in Fig. 13. A parametric form of the Birch-Murnaghan equation of state was used to fit the compressibility along each of the unit-cell axes (for a 3 , b 3 and c 3 ). The density as a function of pressure was fitted using the same equation of state used to fit P(V) points. The fitted parameters are: 0 = 1452.7 (8) kg m À3 , K 0 = 24.2 (3) GPa and K 0 0 = 4.11 (8).
In -ammonium carbamate, the a and b axes compress at a much slower rate than the c axis (Fig. 13). The compression seen along the a and b axes is caused by the similar hydrogenbond geometry in both directions from the R 2 2 (8) ring motifs of the carbamate ions. The structure is much less stiff perpendicular to the planar layers, in which the planar layers of carbamate ions are linked by ammonium tetrahedra halfway between. This can be seen in the compression of the mediumrange hydrogen bonds that connect the layers of carbamate ions via ammonium tetrahedra; the O1Á Á ÁH1 and O2Á Á ÁH2 hydrogen bonds compress by 5.32 and 6.65%, respectively, over a pressure range of 0 to 10.025 GPa. View of the DFT-simulated structure of -ammonium carbamate parallel to the b axis. The grey dashed lines are HÁ Á ÁO hydrogen bonds and the solid black line is the unit cell.

Figure 12
Comparison of the hydrogen-bond motif of -ammonium carbamate in the ab plane from (a) an X-ray diffraction experiment (Kuhn et al., 2007) and (b) the zero-pressure structure from CASTEP with symmetry operators overlain. The relaxed structure shows the same hexagonal-ring motif present in the -phase.

Table 7
A comparison of the unit-cell parameters of -ammonium carbamate from the DFT athermal zero-pressure relaxation using CASTEP and a powder X-ray diffraction (PXRD) experiment at room temperature (Kuhn et al., 2007).
Formula (NH 4   Fits to equation (S5) were achieved using only those pressure values that were greater than 0 GPa. For the axial values, the cubes of the axes were fitted.

Thermodynamics
Since the total internal energies of both structures have been simulated using DFT as a function of pressure, it is possible to assess the relative stability of both phases in order to determine where a phase change may occur as a function of pressure.
The relative stability of the two phases is assessed by calculating the enthalpy (H) as a function of pressure for each phase, calculated using H = U + PV; since the simulations are athermal the enthalpy is equal to the Gibbs free energy. The enthalpies are compared in Fig. 14. The -polymorph has the lowest enthalpy of the two phases and is hence the stable phase at zero-pressure in the athermal limit, with a difference of 399.33 J mol À1 with respect to the -polymorph. When increasing the pressure, the enthalpies of the two phases cross; the transition pressure between the two phases, where ÁH = 0, is $0.4 GPa. The -polymorph has the lowest enthalpy above $0.4 GPa and is then the stable phase in the athermal limit. This trend, of the -polymorph having an increasing enthalpy difference relative to the -polymorph, continues up to pressures of $2 GPa. However, as the -phase starts to undergo subtle structural changes above this pressure, the trend reverses and the enthalpy difference between the two phases begins to decrease.

Thermoelastic behaviour
The presence of negative linear compression in -ammonium carbamate is rather uncommon, since most materials compress in all directions when hydrostatically stressed. It seems rather curious that the negative linear compression behaviour calculated to exist in the -phase is not seen in -ammonium carbamate, since both phases consist of the same centrosymmetric Calculated enthalpy difference between -ammonium carbamate (normalized to zero) and -ammonium carbamate (red squares). Enthalpies were output from the dispersion-corrected simulations of CASTEP, with the lines a simple polynomial fit (solid) and a linear fit (dashed) to the points.

Figure 15
Projections of the thermal expansion coefficient representation surface in the yz (top), xz (middle) and xy (bottom) planes at 60, 120 and 180 K. Solid green lines are positive values. Corresponding projections of the structure of -ammonium carbamate are shown on the right.
-structure, a 'wine-rack' motif that forms between the carbamate ions on the mirror planes and the ammonium ions that form halfway between these layers would cause the carbamate ions on the mirror planes to squeeze closer together, with the ammonium ions moving further apart to accommodate this change. Further experimental work, in particular a high-pressure study, is needed to verify the existence of the calculated NLC seen in -ammonium carbamate, whether a phase transition would occur before the NLC region and to investigate the phase relationship between the two known phases. The measured thermal expansion correlates well with the calculated high-pressure behaviour of -ammonium carbamate. The high degree of anisotropy under pressure in both these polymorphs is due to details in the structure: -ammonium carbamate's b axis is the most compressible, due to few hydrogen bonds forming between layers of ammonium tetrahedra and carbamate ions, which run parallel to the c axis, and hence provide little resistance to compression; -ammonium carbamate is most compressible along c due to planar layers of carbamate ions running in chains in the ab plane, with ammonium tetrahedra forming the interlayer halfway between these parallel chains. The a and b axes have lower but similar compressibility due these to planar layers of carbamate ions.
In terms of thermal expansion, the direction of greatest thermal expansion is perpendicular to the corrugated sheets of alternating ammonium tetrahedra and centrosymmetric carbamate ionic dimers, whereas the minimum thermal expansion is in the plane of these alternating sheets. The intermediate value is parallel to stacks of centrosymmetric carbamate ions running along c, connected by chains of ammonium tetrahedra. The limited bonding between the corrugated sheets running along the b axis provides little resistance to expansion, with a three-dimensional bonding network running in the planes of these sheets resisting expansion (Fig. 15). Based upon the DFT calculations, we would expect that -ammonium carbamate should have a smaller thermal expansivity than -ammonium carbamate, due to its lower compressibility.
When the two phases of ammonium carbamate are compared against thermal expansion measurements of other compounds in this ternary system, the high degree of anisotropy is comparable. The ammonium carbonates are highly anisotropic, typically with one unique axis, and two axes that behave similarly, due to planar layers of interconnected ammonium tetrahedra. Indeed, in ammonium carbonate monohydrate, the c axis has been shown from variabletemperature neutron single-crystal data to exhibit negative linear expansion, whilst the a and b axes exhibit large positive expansivities of similar magnitude (Fortes et al., 2014).

Planetary modelling
The work presented here provides some important values for use in planetary modelling: the density as a function of pressure has been calculated for both structures of ammonium carbamate, accurate crystal structures of both phases of ammonium carbamate for further ab initio simulations have been obtained, and the thermal expansion coefficients ofammonium carbamate are now also known. However, further work is still required to understand the properties of ammonium carbamate in planetary environments. In particular, vibrational spectra are required for unambiguous identification on planetary surfaces, as is a study of the P-T phase diagram under conditions applicable to icy bodies, i.e. high pressure and low temperature.

Conclusion
We report dispersion-corrected density functional theory simulations of bothand -ammonium carbamate, the first neutron powder diffraction study of -ammonium carbamate, together with the first structure refinement at non-ambient temperature, measured at 4.2 K. We also carried out the first thermal expansion measurements of -ammonium carbamate (from 4.2 to 180 K) and have shown how the magnitudes of the axial expansivities are related to the crystal structure. We find evidence for the unusual property of NLC in the a axis of -ammonium carbamate above 5 GPa, with the largest degree of negative expansivity occurring between 9 and 10.5 GPa. This work has provided important parameters to be incorporated into structural models of icy bodies in the Solar System. Further experiments at high pressure are needed to investigate the calculated negative linear compressibility in -ammonium carbamate and the possible phase transition between the two known polymorphs.